Is common practice in Geospatial data science to work with datasets collected at different epoch and or referenced to different reference systems. In this context, the transformation parameters needed to convert data accurately into a more up-to-date reference system are often missed or, when available, are valid for wide areas affecting the accuracy of the transformation. In this briefing notes a simplified approach to derive datum transformation parameters is introduced.
A Datum transformations can be defined as a geometric transformations between two three-dimensional coordinate reference systems. A common method consist in apply a linear transformation in the three dimensional space (x,y,z). A general linear transformation of a vector $x$ to another vector $y$ takes the form
$$y=Mx+t \quad (1)$$Each element of the $y$ vector is a combination of the element of $x$ plus a translation or shift represented by an element of the $t$ vector. The matrix M is called transformation matrix and t is called translation vector. With M being square nonsingular, the inverse relation exist (eq. 2)
$$x = M^{-1}(y-t) \quad (2)$$in which case is called affine transformation.
Limited to the two- and three- dimensional space, six elementary transformations are identified, each representing a single effect. They are geometrically described as: Translation, Uniform scale, Rotation, Reflection, Stretch (Nonuniform scale factors) and Skew (Shear).
<img src="../images/linear-transformation.svg", width="80%">
<img src="../images/helmert.svg", width="80%">
The Helmert 7-parameter transformation, which is an affine (distortion-free) transformation in three dimensions, is extensively used in geodesy to perform Datum transformations. It is applied to geocentric coordinates and can be factored into seven elementary transformations: one uniform scale change, three translations, and three rotations.
Considering a generic point $P$ represented in two orthogonal three dimensional Cartesian spatial reference frames $D_1 (x,y,z)$ and $D_2 (x',y',z')$. For small rotation the direct transformation $P_{D_1} \to P_{D_2}$, is given by
$$ {\begin{pmatrix} x_p \\ y_p \\ z_p \end{pmatrix}}_{{D}_2} = \begin{pmatrix} x^{\prime}_0 \\ y^{\prime}_0 \\ z^{\prime}_0 \end{pmatrix} + (1+k) \begin{pmatrix} 1 & R_z & -R_y \\ -R_z & 1 & R_x \\ R_y & -R_x & 1 \end{pmatrix} \begin{pmatrix} x^{\prime}_p \\ y^{\prime}_p \\ z^{\prime}_p \end{pmatrix}_{{D}_1} \quad (3) $$It is common in Geodesy to separate altimetry from planimetry, in this scenario the affine transformation (eq. 3) can be simplified to a plane roto-translation with isotropic scale variation which requires only four parameters: one scale factor $(\lambda)$, one rotation $(\alpha)$, two translations $(x'_0, y'_0)$. The direct (eq. 4) and inverse (eq. 5) transformations are expressed by:
<img src="../images/plane.svg", width="80%">
To estimate the four parameters $(\lambda, \alpha, x_0, y_0)$ at least two planimetric coordinates in the two systems are needed. However, if more positions are available a least square method (Fitting) can be used solving the linear system:
$$ \left\{ \begin{array}{l l} x'_0 + a x'_p + b y'_p - x_p = 0 \\ y'_0 + a y'_p - b x'_p - y_p = 0 \end{array} \right. \quad (6) $$with:
$$ a = \lambda \cos \alpha $$$$ b = \lambda \sin \alpha $$The linear system (6) can be solved knowing at least two points in $(D_1,D_2)$ once estimated the four unknown parameters $(a,b,{x'}_0,{y'}_0)$ it is possible to derive the rotation angle $\alpha$ and the scale factor $\lambda$ by:
$$ \left\{ \begin{array}{l l} \lambda = \sqrt{a^2 + b^2} \\ \alpha = \arctan \frac{b}{a} \end{array} \right. \quad (7) $$The relation expressed in (eq. 6) can be used in two different ways:
Firs we need to generate a proper test dataset, to do so we'll use a combination of pyproj and numpy. Staring from the data used in the Working with coordinates - Datum-transformation example, we generate a series of 50 random points in two different DATUM :
In [ ]:
#import the pyproj and numpy library
import pyproj
import numpy as np
# set a reference point P with coordinates:
P = (-70.93931369842528, 43.13567095719326)
# define projection UTM 19 N:
# UTM zone 19, WGS84 ellipse, WGS84 datum, defined by epsg code 32619
p1 = pyproj.Proj(init='epsg:32619')
#Find UTM coordinates for the point P(-70.93931369842528,43.13567095719326)
x1, y1 = p1(P[0],P[1])
# define projection: UTM zone 19, Clarke 1866, NAD27 datum
p3 = pyproj.Proj(init='epsg:26719')
# transform the UTM coordinates for the point P to projection 3 coordinates.
x3, y3 = pyproj.transform(p1,p3,x1,y1)
# generate a set of random points in the range of 100 meters from P1
# note: we use a fake altitude to perform a 3D transformation
# the value of 6371 is the ray of the spheroid in km
xrand = (np.random.random_sample((50,))*100)+x1
yrand = (np.random.random_sample((50,))*100)+y1
zrand = (np.random.random_sample((50,))*10)+(6371*1000)
xrand,yrand,zrand
# transform the UTM coordinates for the points [xrand, yrand] to the projection 3 coordinates.
x, y = pyproj.transform(p1,p3,xrand[:],yrand[:])
In [ ]:
# now generate 2 dataframes to store the x,y,z coordinates in the two different DATUM
# and save the reults in a space delimited text file
import pandas as pd
d1 = pd.DataFrame(np.array([xrand,yrand,zrand],dtype=np.float).T, columns=['x','y','z'])
d2 = pd.DataFrame(np.array([x,y,zrand],dtype=np.float).T, columns=['x','y','z'])
d1.to_csv('d1.csv', index=False, header=False, sep=" ")
d2.to_csv('d2.csv', index=False, header=False, sep=" ")
In [ ]:
d1
(transformation performed using pyproj)
In [ ]:
d2
We'll first select the first 10 points in both dataframe and use them to estimate the transformation parameters. Then we'll use the other 40 points in $d1$ as input for the transformation. Finally compare the results with the output of pyproj.
In [ ]:
#d1,d2 subsample
d1s=d1[:11]
d2s=d1[:11]
d=d1[11:]
# save to file
d1s.to_csv('d1s.csv', index=False, header=False, sep=" ")
d2s.to_csv('d2s.csv', index=False, header=False, sep=" ")
d.to_csv('d.csv', index=False, header=False, sep=" ")
In [ ]:
from transform import conforme
In [ ]:
res_conforme = conforme(gcpD1='d1s.csv', gcpD2='d2s.csv', knowD1='d.csv', output='conforme.csv')
res_conforme
In [ ]:
from transform import affine
In [ ]:
res_affine = affine(gcpD1='d1s.csv', gcpD2='d2s.csv', knowD1='d.csv', output='affine.csv')
res_affine
In [ ]:
from transform import helmert
In [ ]:
res_helmert = helmert(gcp1='d1s.csv', gcp2='d2s.csv', inputf='d.csv', output='helmert.txt')
res_helmert
In [ ]:
d2[11:][['x','y']]
In [ ]:
delta_conforme = (d2[11:]['x'].values - res_conforme[:,0], d2[11:]['y'].values - res_conforme[:,1])
delta_conforme
In [ ]:
delta_affine = (d2[11:]['x'].values - res_affine[:,0], d2[11:]['y'].values - res_affine[:,1])
delta_affine
In [ ]:
delta_helmert = (d2[11:]['x'].values - res_helmert[:,0], d2[11:]['y'].values - res_helmert[:,1])
delta_helmert